function [psi,phi] = cmsc2(psi,phi,Q,in,mu,w,dx,g)
% function [psi,phi] = cmsc2(psi,phi,Q,in,mu,w,dx)
%   Inputs:
%           psi     -- fine mesh angular flux
%           phi     -- fine mesh scalar flux
%           Q       -- fine mesh external source
%           in      -- inpute structure
%           mtt     -- material index for each fine mesh
%           mu      -- angle set
%           wt      -- angle weight
%           dx      -- fine mesh delta x
% 
%   This function is a slight variation on CMSC.  Instead of computing the
%   partial currents for a cell and using those in the balance equation
%   (via the discontinuity factors) we instead use the left and right
%   currents across a coarse mesh edge.

% allocate
sigt  = zeros( length(in.xfm), 1 );   % homogenized total cross-section
sigs  = zeros( length(in.xfm), 1 );   % homogenized scatter cross-section
s_c   = zeros( length(in.xfm), 1 );   % volume-averaged external source
phi_c = zeros( length(in.xfm), 1 );   % coarse mesh scalar flux
psi_c = zeros( length(in.xcm), 2 );   % coarse mesh edge angular flux
Jp    = zeros( length(in.xcm), 1 );   % coarse mesh edge rightward (NOT J+ as in description)
Jm    = zeros( length(in.xcm), 1 );   % coarse mesh edge leftward (NOT J- as in description)
Jref  = zeros( length(in.xfm), 2 );   % reference coarse mesh total incoming and outgoing partial current
Jhom  = Jref;                         % homogenized coarse mesh total p cur
fp    = Jp;                           % discontinuity factor for right going
fm    = fp;                           % df for left going
q_c   = zeros( length(in.xfm), 1 );   % coarse mesh emission density
Jp_c = Jp; 
Jm_c = Jm;
% coarse mesh delta x's
h = in.xcm(2:end)-in.xcm(1:end-1);

s = zeros(sum(in.xfm),1);  % the angle-integrated ext src (i.e ausume isotropic)
for i = 1:length(s)
    s(i) = sum( Q(i,:,1)'.*w(:) );
end
Jp(1) = partcur( 1, psi, 1, mu, w, in );
Jm(1) = partcur(-1, psi, 1, mu, w, in );
for i = 1:length(in.xfm)  
    % note, the homogenization is meaningless until I implement a way to
    % specify coarse meshes that contain heterogeneities
    idx1     = 1 + sum( in.xfm(1:(i-1)) );            % lower index
    idx2     = sum( in.xfm(1:(i  )) );                % upper index   
    Jp(i+1)  = partcur( 1, psi, idx2+1, mu, w, in );  % J+ -->
    Jm(i+1)  = partcur(-1, psi, idx2+1, mu, w, in );  % J- <--
   % Jref(i,1)= Jp(i)+Jm(i+1);
   % Jref(i,2)= Jm(i)+Jp(i+1);
    phi_c(i) = sum( dx(idx1:idx2)'.*phi(idx1:idx2,1) ) / h(i); % coarse mesh phi
    sigt(i)  = sum( dx(idx1:idx2)'.*phi(idx1:idx2,1)*in.data( in.mt(i), 1 ));
    sigt(i)  = sigt(i)/(phi_c(i)*h(i)); % total cross-section
    sigs(i)  = sum( dx(idx1:idx2)'.*phi(idx1:idx2,1)*in.data( in.mt(i), 5 ));    
    sigs(i)  = sigs(i)/(phi_c(i)*h(i));
    s_c(i)   = sum( dx(idx1:idx2)'.*s(idx1:idx2,1) ) / h(i);    
end

phi_ref = phi_c;

% compute Qhat

% here, we use an S2-like quadrature.  We want to conserve partial
% currents.  In S2, we have just one angle going right and one going left.
% The partial current J+ = w*mu*psi, or psi = J+/w/mu, and likewise for the
% leftward J-.  We correct Q such that psi(i+1)=J+/w/mu=f(psi(i),Q).

% using an S2-like quadrature
mu_c = 0.00000001;%0.5773502691;
wt_c = 1;

% We have reference partial currents going left and right
% convergence parameters
eps_phi = 1e-6; 
max_it  = 1000;
err_phi = 1;    
it = 0;
% Begin coarse mesh source iterations
NCM = length(phi_c);

ef = zeros(NCM,1);
for i = 1:NCM       % left-to-right
    ef(i)  = exp(-h(i)*sigt(i)/mu_c);
end
while (err_phi > eps_phi && it <= max_it )
    % Save old scalar flux
    phi0 = phi_c; 
    % Update sources
    for i = 1:NCM
        q_c(i) = s_c(i) + sigs(i)*phi_c(i);
    end
    % Perform sweeps
    for i = 1:NCM       % left-to-right
        psi_c(i+1,1) = 0.5*(q_c(i))/sigt(i)*(1-ef(i)) + psi_c(i,1)*ef(i);
    end
    for i = NCM:-1:1  	% left-to-right
        psi_c(i,2)   = 0.5*(q_c(i))/sigt(i)*(1-ef(i)) + psi_c(i+1,2)*ef(i); 
    end    
    Jp_c(1) = psi_c(1,1) * mu_c * wt_c;
    Jm_c(1) = psi_c(1,2) * mu_c * wt_c;
    for i = 1:NCM
        Jp_c(i+1)  = psi_c(i+1,1) * mu_c * wt_c;
        Jm_c(i+1)  = psi_c(i+1,2) * mu_c * wt_c;
    end
%     for i = 1:NCM
%         Jhom(i,1)= Jp(i)+Jm(i+1);
%         Jhom(i,2)= Jm(i)+Jp(i+1);        
%     end
    if ( it == 0 ) 
       fp(2:end) = Jp(2:end)./Jp_c(2:end); % this indexing assumes Jp(1)=0,
       fm(1:end-1) = Jm(1:end-1)./Jm_c(1:end-1); % true for vacuum only
    end
    % Update phi from neutron balance
    for i = 1:NCM
       phi_c(i) = ( q_c(i)*h(i) + ...
                    (fp(i,1)*Jp_c(i,1)+fm(i+1,1)*Jm_c(i+1,1)) - ... % total incident current
                    (fm(i,1)*Jm_c(i,1)+fp(i+1,1)*Jp_c(i+1,1)) ) ... % total  outgoing current
                  / ( sigt(i)*h(i) );
    end
    % Update error and iteration counter
    err_phi =  max(  abs(phi_c-phi0)./phi0 );
    it = it + 1;
end

for i = 1:length(in.xfm)  
    % within cell indices
    idx1 = 1 + sum( in.xfm(1:(i-1)) ); % lower index
    idx2 = sum( in.xfm(1:(i  )) );     % upper index   
    phi(idx1:idx2,1)   = phi(idx1:idx2,1)  * phi_c(i)/phi_ref(i);
end
%figure(2),plot(xc,phi_ref,'k',xc,phi_c,'r--')
%lala=1;
end

function J = partcur( flag, psi, i, mu, w, in )
    if ( flag == 1 )
        idxs = in.ord/2+1:in.ord;
    else
        idxs = 1:in.ord/2;
    end
    J = sum( abs(mu(idxs)) .* w(idxs) .* psi(i,idxs)' );
end

